rm(list=ls())

library(rjags)
load.module("glm")
load.module("lecuyer")
load.module("dic")

jags_data <- "./data/for_jags.Rdata"
jags_model <- "bay_mfa.bug"


# LOAD  
load(jags_data)
cmpemp <- for_jags$cmpemp
phi.select <- for_jags$phi.select 
theta.select <- for_jags$theta.select
rho.select <- for_jags$rho.select


# DATA / CONSTRAINTS
mrg.y <- as.matrix(cmpemp[,c("lr1", "lr2", "lr3", "lr4", "lr5", "lr6", "lr7", "lr8", 
	"lr9", "lr10", "lr11", "lr12", "lr13", "eu1", "eu2", "eu3")] )
mrg.y <- apply(mrg.y, 2, scale, scale=F)

l0 <- c(0,0)
L0 <- diag(2) * 0.01
e0 <- f0 <- 0.001

sigmainv0 <- for_jags$sigmainv0

phi <- matrix(NA,nrow(phi.select),2)
phi <- rbind(phi,for_jags$mu0)

mrg.ce <- as.numeric(cmpemp$ce)

# INITS
inits <- list( 
	list(".RNG.name"="lecuyer::RngStream", ".RNG.seed"=234),
	list(".RNG.name"="lecuyer::RngStream", ".RNG.seed"=104),
	list(".RNG.name"="lecuyer::RngStream", ".RNG.seed"=90),
	list(".RNG.name"="lecuyer::RngStream", ".RNG.seed"=110),
	list(".RNG.name"="lecuyer::RngStream", ".RNG.seed"=919),
	list(".RNG.name"="lecuyer::RngStream", ".RNG.seed"=331) )

# ESTIMATE 
jags <- jags.model(jags_model,  inits = inits, data = list(
								'mrg.y' = mrg.y,
								'N' = nrow(mrg.y), 'L' = ncol(mrg.y),
								'mrg.it' = cmpemp$it,
								'phi.z' = phi.select$z, 
								'phi.w' = phi.select$w,
								'phi.IT' = nrow(phi.select), 
								'phi' = phi, 'mrg.ce' = mrg.ce,
								 'CE' = max(mrg.ce),
								'sigmainv0' = sigmainv0,
								'l0' = l0, 'L0' = L0, 
								'e0' = e0, 'f0' = f0,
								'd' = c(1,1) ), n.chains = 6)

update(jags, 10000)
out_jags <- coda.samples(jags, c('phi', 'lambda', 'sigma', 'K', 'deviance'), 20000, thin=10)

save(out_jags,file="./estimates/out_jags.Rdata")

